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Abstract —To ensure safety in confined environments such as 
mines or subway tunnels, a (wireless) sensor network can be 
deployed to monitor various environmental conditions. One of 
its most important applications is to track personnel, mobile 
equipment and vehicles. However, the state-of-the-art algorithms 
assume that the positions of the sensors are perfectly known, 
which is not necessarily true due to imprecise placement and/or 
dropping of sensors. Therefore, we propose an automatic ap¬ 
proach for simultaneous refinement of sensors’ positions and 
target tracking. We divide the considered area in a finite number 
of cells, define dynamic and measurement models, and apply 
a discrete variant of belief propagation which can efficiently 
solve this high-dimensional problem, and handle all non-Gaussian 
uncertainties expected in this kind of environments. Finally, we 
use ray-tracing simulation to generate an artificial mine-like envi¬ 
ronment and generate synthetic measurement data. According to 
our extensive simulation study, the proposed approach performs 
significantly better than standard Bayesian target tracking and 
localization algorithms, and provides robustness against outliers. 

Index Terms —confined environments, tunnels, sensor network, 
simultaneous localization and tracking, belief propagation, hid¬ 
den Markov model, ray tracing, time of arrival. 


A wireless sensor network (WSN) can be deployed across 
the area to monitor the environmental conditions such as 
stability, temperature and gas levels. The information obtained 
from the sensors can be used to control the ventilation system, 
and determine the unsafe areas and rescue paths. Beyond this 
ability, a WSN can be used to track the personnel, mobile 
equipment and vehicles. The problem is very challenging due 
to the unavailability of GPS signals and the characteristics 
of the propagation environment. The knowledge of the last 
location of an employee is especially important in the af¬ 
termath of accidents such as a wall collapse, explosion, or 
water inundation, but can be also used for task optimization, 
production monitoring and traffic management. For instance, 
according to the MINER act Q, created in response to the 
many mine tragedies in the United States during 2006, the 
emergency response plan “shall provide for above-ground 
personnel to determine the current or immediately pre-accident 
location of all underground personnel”. This problem, that 
also exists in many other conhned environments, is the main 
motivation behind the work reported in this paper. 


I. Introduction 


B. Related Work 


A. Background and Motivation 

A conhned environment represents a constrained and 
irregularly-shaped area, consisting of a series of tunnels or pas¬ 
sages that connect different rooms or halls. Typical examples 
are underground mines, caves, steel factories and subways. In 
these environments, the working conditions may be hazardous 
due to the possibilities of traffic accidents, machine collisions, 
wall collapses, hres and explosions. These environments re¬ 
quire continuous monitoring using sensors deployed all over 
the area. The sensors may be wired and connected to control 
rooms, but to improve the safety and reduce operational costs, 
recently the industry is developing robust wireless communi¬ 
cation systems for this kind of environments 
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Contemporary techniques for localization and tracking in 
confined environments are very basic. They are typically based 
on manual reporting of the employee’s location using paging 
phones or video surveillance &0- Moreover, there are few 
proposals in the literature, based on fingerprinting |j^-| 
trilateration |[^, fT^ , centroid and Bayesian filtering | 

m. 

More specifically, in a hngerprinting technique was 
proposed, in which seven relevant parameters (including mean 
excess delay, total received power, and delay spread) were 
learned offline from wideband impulse responses measured 
at hundreds of locations. Then, these 7D vectors were used 
as the input to an artificial neural network pattern-matching 
algorithm. The measurements were conducted in a gallery of 
the CANMET mine, a former gold mine located in Quebec, 
Canada. This method was then improved in 0 by using more 
receivers with known positions. The fingerprinting techniques 
Gg, in), based on WiEi signals, have been also applied 
in subway tunnels in Seoul, S. Korea. The main problem of 
these algorithms is that they are not well suited for dynamic 
propagation environments (e.g., caused by movement of heavy 
machinery) in which the fingerprints have to be updated very 
frequently. 

In ig, fl^ , ultra-wideband (UWB) measurements were 
used for positioning. They were motivated by a high ranging 
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accuracy in cluttered environments and low-cost implementa¬ 
tion of the communication system. To solve the trilateration 
problem, many types of algorithms have been applied, in¬ 
cluding linearized least-squares, Gauss-Newton and bounding- 
box methods. The measurements were performed in the same 
environment as the one studied in ||^. The main drawbacks of 
these algorithms are that the sensor nodes have to be precisely 
deployed and maintained, and that the algorithms are sensitive 
to outliers. 

In p3j , a centroid algorithm was proposed, in which the 
miner’s location was found by averaging the coordinates of 
the detected anchors. The algorithm is a part of a structure- 
aware self-adapting (SASA) WSN, which is capable of detect¬ 
ing structure variations caused by mine collapses. The main 
problem of this approach is that it requires a high density of 
uniformly deployed sensor nodes. 

In p^ , Bayesian point-mass (grid-based) hltering was 
applied to track mine vehicles. The main goal was monitoring 
and control of ore extraction from the draw points in a mine in 
Australia. Since the draw points are very close to each other, 
high tracking accuracy is required. The main problem of this 
approach is that it requires many grid points in order to obtain 
sufficiently accurate estimates. 

Finally, the results of the measurement campaign |T5), 
carried out in a basement tunnel of Linkoping university and 
an iron-ore mine in Kiruna, Sweden, indicated that UWB time- 
of-arrival (TOA) allows very accurate ranging in line-of-sight 
(LOS) and non-LOS (NLOS) scenarios caused by thin obsta¬ 
cles. However, if the direct path is blocked by a thick tunnel 
wall, the TOA-based ranging leads to a relatively large bias. 
Moreover, the analysis showed that NLOS conditions can¬ 
not be accurately discriminated from LOS conditions, which 
means that (Bayesian) soft-decision algorithms are required 
for accurate ranging and positioning in these environments. 

The previously described state-of-the-art algorithms assume 
that the positions of the sensors are perfectly known, which is 
not necessarily the case due to imprecise placement and/or sen¬ 
sor drops caused by vibrations or wall collapsesj^One possible 
solution to this problem is to manually and periodically verify 
that the sensors positions are correct. However, this approach 
may be too costly and even infeasible in some areas due to 
the on-going activities. 

C. Technical Contributions 

In this paper we propose an automatic approach for target 
tracking with uncertain sensor positions, which involves both 
simultaneous rehnement of the sensors position estimates (lo¬ 
calization) and target tracking (SLAT). Our specihc technical 
contributions are as follows: 

• We divide the considered area into a hnite number of 
cells, and dehne appropriate dynamic and TOA mea¬ 
surement models that take into account the quantization 
effects associated with this division. 

'Although probably not available in confined environments nowadays, we 
also envision that uncertain sensors’ positions can be an outcome of some 
(cooperative) sensor network localization algorithm |16| , GZl 


• We formulate the localization and target tracking problem 
in a Bayesian setting and apply a discrete variant of 
belief propagation (BP). The resulting proposed algorithm 
(referred to as SLAT-BP) can efficiently handle the high 
dimensionality of the problem and the non-Gaussian 
uncertainties. 

• To demonstrate the performance of our SLAT-BP algo¬ 
rithm, we perform an extensive simulation study using 
synthetic impulse responses obtained from ray-tracing 
simulation of a mine-like environment. Our results show 
that SLAT-BP performs signihcantly better than standard 
Bayesian target tracking and localization algorithms, and 
provides robustness against outliers. 


D. Paper organization 

The remainder of this paper is organized as follows. In 
Section we formulate the problem and dehne the dynamic 
and measurement models. In Section III we propose the 
algorithm for simultaneous localization and tracking, based 
on real-time belief propagation. TOA error modelling using 
ray-tracing simulation and performance analysis are provided 
in Section |IV| Finally, conclusions and proposals for future 
work are provided in Section |V] 


H. System model 
A. Problem formulation 

We consider Ns sensors with hxed 3D positions z„ = 
{zn,!, Zn, 2 , n = 1,2,..., Ng, and one target, with 3D 
position xt = at time t, t = 1,2,..., Nt, 

moving through the conhned area. Fig. [T^ illustrates the 
scenario. The sensors are usually placed on the walls or 
the ceiling, but their positions are not perfectly known. A 
moving target periodically emits a signal (including a unique 
identiher^l^that can be detected by a subset of the sensors, with 
a sampling interval T^. The target is equipped with an inertial 
measurement unit (IMU), so it periodically communicates its 
measured velocity. We also assume that there are one or more 
fusion centers (FCs) (e.g., a computer in a control room or a 
target itself), which have available the priors of the sensors’ 
and target’s positions, and periodically collect measurements 
from the sensors and the target. 

A conhned environment is naturally a continuous 3D space, 
but typically it is irregular and it is impossible to analytically 
describe the shape of its borders. On the other hand, using 
an unrestricted 3D continuous space would lead to decreased 
computational efficiency, and more importantly, signihcant 
loss of performance in that position estimates could end up, 
for example, behind a wall. Therefore, we propose to use 
a discrete 3D space, in which the environment is divided 
into a hnite number of cells. The 3D position of the cell 
Ic = ( 4 . 1 ; 4 , 2 j (c = !,■■■, Nc, where is the number 
of cells) is represented, in Cartesian coordinates, by the 
approximation of its geometrical center. It is thereby assumed 

^That means that our algorithm can be also used for multi-target tracking, 
simply by running the same algorithm multiple times. Otherwise, different 
algorithms, e.g., with data association |18| , would be necessary. 
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Fig. 1: Illustration of target tracking in a confined environment: (a) deployment 
of the sensors and a target track, (b) division of the area into 20 cells. 


that the PCs have available a detailed floor plan of the whole 
area. The cell size must be chosen based on a trade-off 
between computational complexity and performance, and it is 
preferable that all the cells have approximately the same size. 
With this model, in which z„ and Xj are discrete variables, 
our goal is to identify in which cells the target and sensors 
are located. This approach also facilitates the application of 
belief propagation (see Section III-B| l without applying Monte 
Carlo or other approximations. Fig. illustrates a confined 
environment divided into cells. 

Finally, we assume that the prior knowledge of the sensors 
and the target positions is defined by probability mass func¬ 
tions (PMFs) Pn,o{^n,o), and po(xo)- While the target’s prior 
represents information about its initial position, the priors of 
the sensors’ positions represent the available information on 
the real positions of the sensors. In practice, this information 
may be provided by the deployment team, and may be 
imperfect due to imprecise deployment or drops of sensors. 


B. Dynamic Model based on IMUs 

We assume that the target is equipped by an IMU, consisting 
of a 3-axis accelerometer, and a 3-axis gyroscope, that aids 
the tracking. IMUs are relatively insensitive to environmental 
conditions, and a good choice in the motivating applications. 
The IMU periodically provides (each Tg second) a velocity es¬ 
timate in all 3 directions Vj = (um, ft,a)- The coordinate 

in the K-th dimension (k = 1,2,3) of the measured velocity 
at time t can be written as a function of the target’s positions 
as: 


Vt,K = 




Tg 






( 1 ) 


^ PqW = Unif(w«; 

2 s ^ s 


where ^ and are samples (at time t, and dimension 
k) of the quantization and measurement noise, drawn from 
the uniform Unif(-), and the Gaussian distribution N{ ), 
respectively. The parameter represents the variance of the 
measurement noise of the IMU, and —D/Tg, D/Tg represent 
the boundaries of the uniform distribution. The quantization 
noise is introduced into the model here because we have to 
relate the continuous quantity u* „ with the discrete quantity 
— Since the boundaries of the irregularly-shaped 

cell are difficult to compute, we approximated the cell with 
a circumscribed cube with length D (D = maxZ?K,c, where 

C,K ’ 

is the maximum possible distance over dimension k in 
cell c). Note that as the cell size D decreases, the quantization 
noise will eventually vanish. 

The distribution of the total IMU noise u = + u™ is 

given by following convolution: 


Puiu) = J Pmiu - u’^)Pq{u‘^)du‘^ (2) 

,, D D 

CX —) - $(u- —), 

where $(m) = 0.5 (l-f erf(M/(CT„v/2)) is the cumulative 
Gaussian distribution. Now, we can model the dynamic^ of 
the target as: 


p(vt,Xt|xt_i) =p(vt|xt,xt_i)p(xt|xt_i) (3) 

cxp(vt|xt,X(_i) 

= H p{vt,K\xt,,i,Xt-i,K) 
k^1,2,3 

- Pui^Xt^l^ {Xt^K Xt—l^K^/Tg)^ 

k=1,2,3 

where we assumed independence between the coordinates, and 
that the target’s mobility model is unknown (i.e., p(xt|xt_i) cx 
l)j^ Recall also that Xj is discrete, so Q gives us information 
about the cell of the target at time t, given the cell of the target 
at time t — 1, and the measured velocity v^. 

For our framework (see Section |III-A| ), it is convenient to 
formally define the “dynamics” of the sensors: 




(4) 


where 6{-) is the Dirac delta impulse, which enforces that the 
sensors are static (Zn,t = Zn,t-i = z„). 

Finally, we note that, if po(xo) is very informative (e.g., the 
target’s initial cell is given), the target’s dynamic model (given 
by ([^) already constitutes sufficient information for a tracking 
algorithm known as dead reckoning |19|, but this approach 
would suffer from error accumulation over time. Therefore, we 
need to use a WSN which will provide periodic measurements 


^In some of the literature, the quantity p(vt,xt|xt_i) is denoted hy 
p(xt|xt_i). We prefer to write out the measured velocity explicitly. 

'^Even without this assumption, the IMU typically provides much more pre¬ 
cise information than the mobility model, so p(vt|xt, Xi_i)p(xt|x(_i) ^ 
p(vt|xt,Xt_i). 
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w.r.t. their positions. 


C. TOA Measurement Model 

We assume that round-trip TOA (RT-TOA) measurements 
are obtained, at each time slot, by a subset of the sensors in 
proximity to the target. We decided to use RT-TOA (instead 
of one-way TOA), in order to avoid the need for clock 
synchronization between the sensors and the target. We did 
not consider received signal strength (RSS) since the distance 
estimates will be highly erroneous, due to the severe multipath 
in confined environments, as shown in fig, 1^, 1^. RT- 
TOA can be obtained using many techniques p2) , but we 
assume that a simple thresholding is performed. More exactly, 
the RT-TOA is taken to be the arrival time of the first multipath 
component in the measured impulse response that exceeds 
a predefined threshold. Note that a signal should have very 
large bandwidth (WB or UWB), in order to provide high TOA 
resolution j^ . 

The measured distance between sensor n and the target at 
time t, can be written as; 


dt,n — dpT — 


■ 


■ Wt 


(5) 


“t,n 


. Pq{w^) = Unif(w«; 0, Ds/l), ~ Pmiw^), 


where dpT is a known bias caused by processing time on a 
target, is measured TOA, c = 3 • 10® m/s is the speed 

of light, and wf „, are samples of the quantization noise, 
and measurement noise, respectively. The distribution of the 
quantization noise is not available in a parametric form, so 
we choose the least informative (i.e., a uniform) distribution 
to keep the algorithm tractable. 

While the measurement noise depends on many factors 
(such as thermal noise, bandwidth, and the quality of the 
sensors), the bias arising from multipath propagation in NLOS 
conditions is usually the most critical source of error. The most 
common approach p4) , | [25j is to identify NLOS measure¬ 
ments, and discard them, or alternatively mitigate the effect 
of the multipath bias. However, these techniques require an 
identification of NLOS conditions, which cannot be always 
accurately done (see, for example, fB), @). Therefore, we 
prefer a model that does not require NLOS identification, but 
only knowledge of the probability of having NLOS. 

According to previous results Q, we can roughly model 
line-of-sight (LOS) measurement noise with a Gaussian dis¬ 
tribution, and NLOS noise (in which the walls block the 
direct path) with a Weibull or an exponential distribution. 
However, in severe multi-path environments, it is expected 
that NLOS noise has multiple modes (see also Section IV-A| l. 
Therefore, we use a Gaussian mixture (GM), which is capable 
of approximating arbitrary probability distributions. Moreover, 
we also need to take into account other sources of NLOS, 
such as vehicles and machinery. Since these objects are usually 
dynamic, made of different materials, and have different sizes 
and thicknesses, it is difficult to model their effect. Therefore, 
we assume that this error is uniformly distributed (14), but 
that it appears with very small probability. In total, the model 


for Pw{w'^) is then given by the following mixture; 

= Flos ■ A/'('u;"';0,cr^ o) + (6) 

FnloS * ^ ^ ^w^i) 

PoBS-Unif(u;”^;0,D„,ax), 


where CTu, o is the standard deviation of the LOS component of 
the noise; pw^i, o'w,i are the weights, means and standard 
deviations of the NLOS noise caused by tunnel walls; Nm is 
the number of GM components; and Omax is the maximum 
distance error. Plos^ Pnlos and Pqbs are the probabilities 
of LOS, NLOS caused by tunnel walls, and NLOS caused by 
other obstacles, respectively (with Plos = 1—Pnlos~Pobs)- 
While Pnlos and Pqbs can be approximately estimated by 
examining the floor plan of the deployment area, the GM 
parameters Pw,i and aw,i) can be estimated by applying 

the expectation-maximization (EM), generalized EM, or the k- 
means algorithm Chapter 9] on training samples. 

,,9 


The distribution of the total noise w = w‘^‘ 
by the following convolution; 


+ IS given 


Pw{w) = J Pm{w - 'w‘^)Pq{w‘‘)dw‘‘ (7) 

= ^ (4>o(u;) - - Pv^)) + 

g{w)PoBS, 


where is the shorthand notation for $(•; i), and 

the distribution g{w) is found by convolution of two uniform 
distributions; 


9{w) 


W / (PmaxPV^), 0 < W < Pa/3 

1/Pniax ; Ds/S < tU < Pmax 

(Pmax ~ tf)/(PmaxP'\/3)j Pmax < W < Pmax 
0, otherwise 


( 8 ) 

where Pj„ax = Pv^ + Pmax- Einally, the likelihood function 
is given by; 


Pidn,t\^t,2n,t) =Pw{dn,t “ \\^t “ Zn.tll)' (9) 


The likelihood functions can be computed for all sensors 
which detect the target. However, we assume that the sensors 
can perform measurements with the target if and only if dn,t < 
c/th where c/th is a predefined sensing radius. The set of 
sensors that perform measurements at time t is denoted by Gf. 
The sensing radius is chosen so as to ensure that a sufficient 
number of sensors can detect the target, but should be small 
enough that the model in (|^ remains valid. 


111. Simultaneous Sensor Localization and Target 
Tracking 

Our goal is to obtain the posterior marginal PMEs (referred 
to as the beliefs), p(xt|ei:t) and p(z„_t|ei:t), of the following 
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Fig. 2: An example of a MRF, with three sensors (n, p, and r) and a target 
at three time instants (t — 1, t, and t + 1). The messages are represented with 
an'ows. 


Note that these potentials, as well as all messages and beliefs 
(defined in further text), are not necessarily normalized. 

B. SLAT via Real-Time Belief Propagation (BP) 

We adapt the standard BP (see ||^ eqs. (8)-(9)]) for our 
real-time and discrete problem. To ensure real-time execution, 
we do not send the messages backward in time (see also Fig. 
1 ^, and since the variables are discrete, we replace integra¬ 
tion with summation. The beliefs (denoted with M„ t(z„ t) 
and Mt(xt), respectively) are initialized with M„o(zn,o) = 
'0n,o(zn,o) Mo(xo) = ^o(xo)- The algorithm (for time 
slot t, t > 0) is summarized below; 

Step 1. Compute the sensor-to-target and target-to-target 
messages: 

W(n.t)->t(xt) = (16) 

TOt_l_>t(xt) = ^ V’i-l.t(xt-l,Xt)Mt_l(X4_l). (17) 

Xt_l 


joint distribution; 

p(xo,...,Xt,Zi, 0 ,---,ZAf,,t|ei:t) oc (10) 

Po(xo) n Pn.o(Zn,o) H jx*/, Z„,t/) • 

n—l...Ns = 

n p(vt/,Xt.|xt/_i) n p{Zn,t'\Zn,t'-l), 

where ei-t is all available evidence up to time t (i.e., measured 
TOAs and velocities). The previous factorization is obtained 
using Bayes’ rule and standard assumptions | |27| , such as 
independence of the measurements/priors and memoryless 
movement. Since the marginalization of ( [T0) i is intractable, 
we resort to message-passing on a graphical model. 


A. Graphical model 

We use an undirected graphical model p 8 | , also known as 
a Markov random field (MRF)j^to represent the factorization 
in ( [Tol l. In ^ MRF, each vertex represents a random variable 
with an associated single-node potential (a local evidence), 
and each edge represents a pairwise potential (a likelihood 
function). An example is shown in Fig. [^ Using the models 
defined in Sections ITB and |II-C| the potentials are given by; 


V’i(xt) = <j 

f Po(xo), if < = 0 
[ 1 , otherwise ’ 

( 11 ) 

V’n,t(Zn,i) — j 

f Pn, 0 l 
1 

<^ 71 , 0)7 if f — 0 
otherwise ’ 

( 12 ) 


’n,i) — 

p(d„,t|xt,z„,t), 

(13) 

V't-i,t(xt_i 

,Xt) = 

:p(vt,Xt|Xt_i), 

(14) 

^n,i —l),(n,i) (^n, 

t—1? 

fc,t) — —l)- 

(15) 


^An equivalent algorithm can be derived using the forward phase of the 
forward-backward algorithm (also known as BCJR) in a hidden Markov model 
(29| . However, we prefer to use a much more flexible framework, valid for 
discrete, continuous and mixed variables. 


Step 2. Update the target’s belief; 

m(„,t)_>t(xt). (18) 

neGf 

Step 3. Compute the target-to-sensor and sensor-to-sensor 
messages; 

/ V _ / r \ 

/ \? (19) 




( 20 ) 


Step 4. Update the beliefs of the sensors: 




Mn,t-l{2n,t)mt^(^ri,t){Zn,t), if ^ € Gf 
otherwise 


( 21 ) 

Step 5. (optional) Compute the estimates using the k- 
nearest neighbour (kNN) approach 


— 




( 22 ) 


x; xtMt{xt) 
E Mtixt) 


(23) 


where ^ and G^^ are the set of k cells with highest beliefs, 
t(z„ () and Mt{Kt), respectively. The special cases k = 
1 and k = Nc, correspond to MAP and MMSE estimates, 
respectively. Note that this phase is optional, since the main 
output of this algorithm (beliefs) have been already computed 
in steps 2 and 4. 

For comparison purposes, we also consider two specific 
instances of the proposed SLAT algorithm; i) Bayesian point- 
mass target tracking (by excluding all target-to-sensor mes¬ 
sages, i.e., mi_j.(„^t)(z„ t) = 1), and ii) Bayesian point- 
mass target localization (by excluding all target-to-sensor 
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and target-to-target messages, i.e., t) (z„_() = 1 and 

= 1). The former one uses the sensors’ priors to 
track the target, while the latter one uses the sensors’ priors to 
locate the target independently in each time slot (i.e., without 
a dynamic model). Note that the target tracking algorithm in 
03 can be considered as a special case of SLAT-BP (although 
their measurement and dynamic models are different). 


C. Implementation Issues 

In this section, we discuss some important issues that can 
arise during the implementation of the proposed SLAT-BP 
algorithm. 

• Complexity of the algorithm: The complexity of the 

SLAT-BP algorithm at time t is 0{\G^ \ N^), since the 
message computations dominates the cost. This complex¬ 
ity is significantly less than that of naive marginalization 
of (101 which would require operations at 

time instant t. Although this is a significant reduction, 
the complexity is still high if there are many cells. The 
complexity can be further reduced by considering only 
beliefs with a probability larger than a predefined belief 
threshold cm Gm < !)■ For example, in ( [Th] ), only cells 
c, which satisfy the following constraint: 

^ >eM/iV„ (24) 

c' = l...Afc 

should be considered for summation. An analog con¬ 
straint is then used for and Denoting the 

number of these cells by ^ and , and = 

max(max„gGj AT^J, the complexity at time t is 

reduced to 0{\Gf\NcN^). 

• Non-synchronized measurements: In Section II-B we 
assumed that the IMU is configured to operate at the 
same rate as the sensors (reporting every Tg second). 
In practice, the rate of the IMU may be much higher 
than that of the sensors. The proposed message-passing 
algorithm can be easily adapted to a situation where the 
rates are different. Assuming that the algorithm operates 
at the IMU rate and that we do not know the rate of the 
sensors’ measurements, we just need to do following at 
each time slot: i) if the sensors’ measurements are avail¬ 
able, we run the algorithm in Section |III-B| and ii) if the 
sensors’ measurements are unavailable, simply exclude 
all sensor-to-target messages (i.e., m.(„,t)_>.t(xt) = 1). 
In other words, we run simultaneous localization and 
dead reckoning in all time slots in which only the IMU 
measurements are available. 

• Routing data to the PCs: All collected measurements 
should be routed to the PCs as soon as they become 
available. Although there are many well-known routing 
protocols p2] , we recommend a hybrid system based 
on a leaky feeder system (LFS) and a wireless mesh 
network (WMN), similar to one installed at a coal mine 
in West Virginia, US Q. LFS consists of a coaxial- 
type cable, which emits and receives radio waves (i.e, 
it behaves as a distributed antenna). It has many power 
supplies, and a backup battery in an explosion-proof 


enclosure. Therefore, all data transmitted by the sensors 
or the target, will be available to FCs using a one-hop 
communication link and without any routing protocol. 
Since LFS typically cannot provide coverage all over 
the deployment area, it should be complemented with a 
WMN, which should consist of the subset of WSN not 
in the vicinity of the LFS cables. In the WMN, sensors 
communicate in a multi-hop fashion using an optimal 
path, computed in real-time. If one or a few sensors fails, 
the system simply recomputes the path. Therefore, this 
system is capable of routing the data as long as there is 
one path between a sensor and a FC or LFS cable. 

• Online calibration: Although the measurement model in 
(|^ provides robustness against dynamic obstacles, it will 
not be good enough if a permanent change is made to 
the environment (e.g., new pillars are formed, or the 
tunnels are extended). In that case, it would be necessary 
to repeat the calibration (especially, to re-estimate the 
GM parameters for the NLOS error model), which is a 
cumbersome task. An alternative, and preferable option, 
is to update the measurement model online using already 
deployed sensors. This is feasible since the fusion center 
knows the current estimates of the sensors’ positions, and 
consequently all inter-sensor distances. The additional 
requirement is that the sensors are randomly deployed 
so that they can provide sufficient statistics for parameter 
estimation. This calibration should be done periodically 
(e.g., once per day), or manually triggered once some 
change in the deployment area is reported. 

IV. Numerical results 

In this section, we analyze the accuracy and the robustness 
of the proposed approach using ray-tracing simulation. 

A. TOA error modelling using ray-tracing 

For TOA error modelling, we decided to use REMCOM’s 
Wireless InSite ray-tracing simulator 1^, 0. Wireless 
InSite is a flexible, powerful tool for accurately predicting 
the effects of the environment on the propagation of electro¬ 
magnetic waves. It models the physical characteristics of the 
environment (including the effects of terrain, urban buildings 
and tunnel features), performs the electromagnetic calcula¬ 
tions, and evaluates the signal propagation characteristics. The 
calculations are made by shooting rays from the transmitters, 
and propagating them through the environment, until they 
arrive at the receivers. The rays interact with environmental 
features via: reflections at object faces, diffraction around ob¬ 
jects, and transmission (penetration) through objects. Wireless 
InSite can provide quantities such as electric and magnetic 
field strength, received signal strength, time of arrival, path- 
loss, delay spread, direction of arrival, impulse response and 
power delay profile. 

We designed an artificial mine tunnel using Wireless InSite 
(see Fig. |^, by creating many small pieces, and connecting 
them together. The dimensions of the tunnel are similar to 
a real mine tunnel |j^, and each piece contains expected 
irregularities. We chose concrete as material for the tunnel 






(a) 


Rx routes 6-9 



(b) 

Fig. 3: Illustrations of the designed tunnel in Wireless InSite: (a) 3D illustration, and (b) deployment of Tx/Rx through the tunnel. The dimensions of the 
tunnel are approx. llOm (length) x 5m (width) x 5m (height) (corresponding to example in Fig.[^). 


walls, which corresponds to the areas in the mine reinforced 
to increase stability. Then, we deployed transmitters (Tx) and 
receivers (Rx) routes along the tunnel, as shown in Fig. 3b 
Our goal is to obtain TOA at each receiver, which has NLOS 
to the transmitter caused by the tunnel walls. LOS links are 
not considered since they provide ground-true estimates (the 
TOA estimates from Wireless InSite do not include other 
sources of errors, such as thermal noise, or limited bandwidth). 
Finally, we set the parameter^ as shown in Table [l] We 
chose the Full-3D Shooting and Bouncing Ray (SBR) method, 
which includes the effect of reflections, transmissions, and 
diffractions on the electric held in 3D environment, without 
any restriction on the object shapes (more details about this 


method, and possible alternatives can be found in 134 Chapter 


15]). Regarding maximum number of reflections, transmis¬ 
sions, and diffractions, taking into account the analyses in lH), 
we chose the values which provide a good trade-off between 
the performance and complexity (i.e., any further increase of 
these values would lead to a negligible difference in the results, 
while the computational time would grow dramatically). We 
used short-dipole antennas with vertical polarization in order 
to ensure a near-omni-directional radiation pattern. Finally, the 
values for input power, transmission line loss, and received 
power threshold were chosen to ensure sufficient range for 
ray propagation. 


“Many other parameters (see I 
values since they are irrelevant 1 


g ), not shown in Table^ are kept to default 
this analysis. 


TABLE I: Main parameters used for Wireless InSite simulations 


Ray-tracing method 

Full-3D SBR 

Antenna (Tx/Rx) 

short dipole 

Polarization (Tx/Rx) 

vertical 

Relative permittivity of the wall 

7 

Conductivity of the wall 

0.015 m-ifl-i 

Waveform (Tx/Rx) 

sinusoidal 

Central frequency 

2.4 GHz 

Input power (Tx) 

15 dBm 

Received Power Threshold (Rx) 

-no dBm 

Transmission line loss (Tx/Rx) 

0 dB 

Altitude of antenna (Tx/Rx) 

1.3 m 

Maximum number of reflections 

8 

Maximum number of transmissions 

1 

Maximum number of diffractions 

1 


From the simulation, we collected 1130 impulse responses, 
and obtained the same number of TOA samples (in the absence 
of noise, the threshold for TOA estimation is set to zero). 
The samples were used to And the distance estimates, and the 
samples of the ranging error (the difference between the true 
and the estimated distance), in order to create a GM model. 
Observing the histogram in Fig. we can see that the error 
is positive (as expected in the absence of the other noise), and 
that it can be well approximated with 5 modes. Moreover, we 
can see that the error is not exponentially distributed, which 















0.045 


TABLE II: Summary of parameters used in simulations 



Fig. 4: Histogram of the NLOS ranging error, and coiTesponding Gaussian 
mixture model. 


Parameter 

Value 

number of cells (A),) 

44 

number of sensors (Ng) 

25 

number of time slots (Nt) 

40 

sampling interval (Tg) 

1 s 

std. deviation of sensors’ positions (as) 

6 m 

prob. of NLOS (tunnel wall) (Pnlos) 

0.17 

prob. of NLOS (obstacle) (Pqbs) 

0.03 

std. deviation of IMU noise (cr„) 

0.5 m/s 

std. deviation of LOS noise (a^jfi) 

1 m 

sensing radius (c/th) 

30 m 

maximum distance error (Umax) 

c/th 

quantization noise (D) 

5 m 

kNN parameter (k) 

2 

belief threshold (cm) 

0.05 

number of Monte Carlo runs (A^mc) 

100 


stands in contrast to results | |35l , for other environments. 
Therefore, we set the number of GM components to Nm = 5, 
and run the k-means algorithm to obtain its parameters]^ 

According to Fig. this model provides a good (but not 
perfect) approximation of the real PDF. We did not include 
other error sources in the NLOS model, since they typically 
are negligible compared to the multipath error (recall that the 
signal-to-noise-ratio is very high for short-range sensing, and 
that we target WB/UWB signals). For the LOS error model, 
we did not analyze the precise impact of the transmit power, 
bandwidth and other sources of the error. We simply assumed 
that it can be well approximated with a zero-mean Gaussian 
distribution with a sufficiently large standard deviation = 
c ■ 3.33 ns = 1 m. 

Finally, we considered a scenario with the same deployment 
as in Fig. but with different types of the obstacles (objects 
created of metal, water, sand, etc.) placed in front of the 
transmitters. We found that the bias caused by these obstacles 
is very uncertain, so our assumption that it is uniformly dis¬ 
tributed (as a least informative option) is reasonable. Assuming 
that these obstacles will block the LOS path in 3% of the 
cases, we increased the size of the NLOS database to 1164 
samples by adding (randomly chosen) samples from the NLOS 
scenarios with obstacles. This database, along with the GM 
model, will be used in the next sections as an input to the 
performance analysis simulations. 

B. Simulation Setup 

We considered the tunnel in Fig. divided in N^. = 44 
cells, and with Ng = 25 sensors randomly deployed in 
these cells. The sensors’ priors were given by Pn,o(^n,o) = 
A/'(z„ o; Iri) ^s) where 1„ (n = l,...,Ns) is the reported 
location provided by the deployment team, and Es = 
diag(CT|, ct|, cr|) {as = 6 m) is the empirical measure of 
the precision of the sensor placement. We assume that we 

^The obtained parameters are not provided since they are valid only for 
this, artificially generated, tunnel. For other environments, the error modelling 
should be repeated. 


know the initial cell of the target, i.e., po(xo) = — li)- 

There are = 40 time slots, and the sampling interval is 
Tg = 1 s. The target position at time t is generated using the 
following mobility model; 

f ht+r^, for 1 <t < Nt/2 + 1, 

* \ hiN^-i-t)+r ,, for Nt/2 + 2 < f < Wt ’ ^ 

where p is a random integer between —1 and 1 (i.e., ry ^ 

Unif{— 1,0,1}), which adds uncertainty to the mobility of the 
target. This model assumes that the target is moving forward 
through the tunnel during the first half of the time, and then 
going backward until the end of the period. Recall that this 
model is not known to our algorithm, but the target’s velocity 
is measured by the IMU. 

We set the remaining parameters for the measurement noise 
as follows: Pnlos = 0.17, Pqbs = 0.03, and (t„ = 0.5 
m/s. The measurements are generated using these models, 
except in the NLOS case, in which we pick a randomly- 
chosen sample from the ray-tracing database. To take into 
account the cell size, we set the quantization noise parameter 
to P = 5 m. Moreover, the sensing radius was set to c/th = 30 
m, the maximum distance error to Umax = dTH = 30 m, 
the parameter for kNN estimation to fc = 2, and the belief 
threshold to cm = 0.05. All results were averaged over at 
least Nmc = 100 Monte Carlo runs. In each run, we generated 
different observations, target tracks and sensor positions. 

All parameters, summarized in Table were used unless 
otherwise stated in the following text. 


C. Performance Analysis 

Our goal is to analyze the accuracy of the proposed SLAT 
and compare it with corresponding tracking and localization 
algorithms (defined in Section III-B| i. The target’s (sensor’s) 
position error is defined as the Euclidean distance between 
the position of the true and the estimated cell of the target 
(sensor). 

We first analyze the root-mean-square error (RMSE) as 
a function of time, as shown in Eig. Regarding the 
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Fig. 5: RMSE as a function of time for (a) the target’s position error, and (b) the sensors’ position error. For the sensors, there is an additional averaging over 
all sensors’ position en'ors. 
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Fig. 6: CDF of (a) the target’s position error, and (b) the sensors’ position error. 


estimates of the target positions (Fig. |^, SLAT provides 
the best performance, followed by tracking, and localization 
which provides the worst performance. The difference between 
tracking and localization is caused by information from the 
IMU. The SLAT performs better than tracking because the 
target uses the information from improved sensors’ PMFs, 
while only the information from initial sensors’ PMFs in 
the tracking algorithm. We also note that the performance of 
all algorithms is the worst when the target is close to the 
tunnel edges (t = 1,21,22,40). This behavior is caused by 
the small number of the sensors in proximity of the target. 
Regarding the estimates of the sensors’ positions (Fig. [5b)l, the 
SLAT consistently improves their estimates comparing with 
localization/tracking which do not update the sensors’ PMFs. 
This can be explained by the fact that the state of the same 
variable (in contrast to tracking) is estimated in each time 


slot. However, the RMSE will not converge to zero because 
the NLOS measurement model is not fully consistent with the 
measurements (see Fig. |^. 

To draw more precise conclusions, we analyze the cumu¬ 
lative distribution function (CDF]j of the target’s and the 
sensors’ position error in Fig. |6 We can see in Fig. 
that the localization, tracking and SLAT algorithms correctly 
estimate the target’s cell in about 42%, 46%, and 53% of 
the considered tests, respectively. Regarding the sensor posi¬ 
tion estimates (Fig. [6b| ), the localization/tracking algorithms 
correctly estimate the sensors’ cells in only 18% of the 
tests, while SLAT does so in about 45% of the tests. More 
importantly, the 95th percentile error of the SLAT (for both 
the sensors’ and target’s RMSE) is at most the half of the 

*The CDF contains discrete steps due to the finite number of possible error 
values. 
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Fig. 7: RMSE as a function of: (a), (d) sensing radius, (b), (e) number of sensors, and (c), (f) standard deviation of the sensors position. Plots with RMSE 
of the target’s position error are shown in the first row, while plots with RMSE of the sensors’ positions error are shown in the second row. 


corresponding localization/tracking algorithms. Therefore, we 
can conclude that SLAT outperforms tracking/localization in 
terms of both the target’s and the sensors’ error. In both cases, 
SLAT is consistently (at any percentile), better than the other 
two algorithms. It is worth noting that these conclusions are 
consistent with our previous results ||T| based on models from 
the CANMET mine. 

We now analyze the effect of the different parameters on 
the performance. The results are shown in Fig. We make 
the following observations; 

• Ejfect of varying the sensing radius (Figs. and 
Increasing the sensing radius consistently leads to lower 
RMSE of the target position estimates (for any of the con¬ 
sidered algorithms), and of the sensor position estimates 
(only for SLAT). This behavior is expected since a higher 
radius allows more sensors to perform measurements, 
which are then used to update both the target’s and 
the sensors’ positions. However, recall that the sensing 
radius can be only increased up to a limit defined by the 
maximum sensing threshold for which the measurement 
model is valid. We also note that the difference between 
the tracking and localization algorithms is decreasing 
with an increasing sensing radius. This can be explained 
by the fact that the IMU provides a smaller proportion 
of the total information when more sensors provide 
measurements. 

• Effect of varying the number of sensors (Figs. an d^): 


Increasing the number of deployed sensors will obviously 
lead to better performance of all considered algorithms. 
Thanks to the measurements from the IMU, the sensitivity 
is not significant for the target position estimates of the 
tracking and SLAT algorithms. Regarding the sensors’ 
SLAT estimates, the sensitivity is extremely low because 
each sensor receives the same number of messages from 
the target regardless of the number of sensors (see Fig.|^. 
Nevertheless, these messages are more informative since 
they are functions of improved target’s beliefs. Finally, we 
can conclude that the SLAT algorithm can be used with a 
lower density of deployed sensors than the corresponding 
localization/tracking algorithms (e.g., if an RMSE of 3 m 
is acceptable, we need 10 sensors in case of SLAT, 15 
sensors in case of tracking, and 22 sensors in case of 
localization). This is especially important if the sensors 
are expensive. 

• Effect of varying the standard deviation of the sensors ’ 
positions (Figs. B Fig. ^): In all analyzed cases, the 
RMSE is nearly a linear function of CTs. However, the 
most important difference between the algorithms is the 
slope of the curve, which is the lowest in case of SLAT 
for both the sensor and the target position estimates. That 
basically means that SLAT is the most useful in scenarios 
in which the sensors are very imprecisely deployed, or if 
most of the sensors drop far from their original positions 
(which may happen in the aftermaths of accidents). 
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Fig. 8: RMSE of the SLAT algorithm as a function of parameters of the contamination noise for (a) the target’s position error, and (b) the sensors’ position 
error. 


Finally, we analyze the robustness of the SLAT algorithm 
against outliers in the distance measurements. Therefore, we 
assume that the measured distance is contaminated with an 
outlier do, with weight Pq- In other words, the contaminated 
measured distance is given by “f where rio is 

contamination noise distributed as; 

Uo(1 - Po)S{no) + PoS{no - do). (26) 

Note that the presence of these outliers is unknown to the 
SLAT algorithm, so it is not part of the model in (|^. These 
outliers may be caused by additional NLOS conditions, inter¬ 
ference, or other sources of errors. According to the results, 
shown in Fig. we hnd the following; i) higher values of 
Po and/or do would lead to increased RMSE, ii) small value 
of Po would not lead to performance degradation regardless 
of the value of do, hi) the sensor position estimates are less 
sensitive to outliers than the target position estimates, and iv) 
for Po > 0.4 and large values of do, the sensor position esti¬ 
mates are worse than the corresponding tracking/localization 
estimates, in which the RMSE is equal to 5.9 m (see Eig. 
even in the presence of outliers. In principle, the algorithm is 
very robust against outliers, especially for small weights Po, 
and this is achieved thanks to the uniformly distributed tail in 
the distribution of the measurement noise (see (|^). The unique 
problem is that SLAT should not be used if these weights are 
large, but this situation should not be expected in reality. If 
the most of the links are contaminated with an outlier (this 
situation can be detected by comparing tracking and SLAT), it 
is necessary to perform re-calibration, as explained in Section 

EE] 


D. Experimental Example 

We consider a small part of the CANMET mine Q 
to test the performance of our algorithm. The tunnel is 
divided into Nc = 14 (non-equal) cells, and in each of 
them there is one sensor (i.e., = No). The sensors’ 


priors are given by p„,o(zn.o) = A/'(z„,o; In, Eg) where 1 „ 
(n = 1 ,..., Ns) is the expected location of the sensors, and 
E 5 = diag(25m^, 25 m^, 0). There are Nt = 30 time slots, 
the sampling interval is Tg = 1 s, and the quantization noise 
parameter is 79 = 6 m. Taking the results from the 
measurement noise follows a Gaussian distribution for LOS, 
and a Weibull distribution for NLOS. More details about the 
considered scenario, and all other parameters, are available in 

Q. 

We analyze the CDE of the target’s and the sensors’ position 
error, shown in Eig. As we can see, SLAT is consistently (at 
any percentile) more accurate than corresponding tracking and 
localization algorithms, which is consistent with our results 
based on ray-tracing (Eig. |^. One minor difference is that 
CDE is smoother than in Eig. 6 , which is caused by variable 
cell sizes. 


E. Other solutions in literature 

Although not for the same type of the environment, there 
are other solutions in literature that try to address a similar 
problem. They are mainly based on hngerprinting techniques 
| 8 K pT], p^ , and SLAT algorithms over continuous space 
I ^-||41[. In general, hngerprinting is able to outperform our 
proposed SLAT method, assuming that enough hngerprints 
are available. However, this approach would require an 
exhaustive calibration that sometimes may be infeasible. An 
additional problem of hngerprinting is that sensors may fall 
from the walls during the training phase (which have to be 
repeated very frequently), making a subset of the hngerprints 
invalid. Regarding other SLAT algorithms, to the best of 
our knowledge, none of them are adapted to non-Gaussian 
measurement models and for conhned areas. Eor example, the 
SLAT methods in |39|-|411 provide the posterior in Gaussian 
form, and use an unrestricted continuous space. Consequently, 
they would either fail in the presence of high levels of NLOS 
bias or provide invalid estimates (e.g., behind the walls). 
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Fig. 9: CDF of (a) the target’s position en'or, and (b) the sensors’ position error. The results are based on models from the CANMET mine. 


V. Conclusions and Future Work 

We presented a novel approach to target tracking in confined 
environments in the presence of uncertain sensors positions. It 
is based on a SLAT-BP algorithm, which can simultaneously, 
and in real-time, refine the sensors’ positions and track the 
target. This algorithm can: i) efficiently solve high-dimensional 
problems, ii) handle all non-Gaussian uncertainties, and iii) 
provide only valid position estimates thanks to the use of cells 
with predefined positions. According to our simulation results, 
both the sensor and the target position estimates are improved 
even after a very short tracking period. Moreover, SLAT-BP 
can be used with a very low density of deployed sensors, and 
can keep performing well even if most of the sensors drop far 
from their original position. Finally, the algorithm is also very 
robust against outliers as long as they appear with reasonable 
probabilities. 

By no means the study in this paper provides the solutions 
for all problems in confined environments. There remain 
many open directions for future work. First, it would be 
useful to provide a method that learns and adapts the ranging 
distribution online. As a result, the SLAT algorithm can be re¬ 
used for a wide variety of environments. Second, a distributed 
implementation of the SLAT algorithm, in which the sensors 
estimate the target’s position by cooperating only with the 
neighboring sensors, may be useful to increase the scalability 
and robustness. Finally, cooperative infrastructure-free self¬ 
localization algorithm would be of high interest for search- 
and-rescue operations. 
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